Multiscale Modelling of Lytic Polysaccharide Monooxygenases

Lytic polysaccharide monooxygenase (LPMO) enzymes have attracted considerable attention due to their ability to enhance polysaccharide depolymerization, making them interesting in respect to production of biofuel from cellulose. The LPMOs are metalloenzymes that contain a mononuclear copper active site that can active dioxygen. However, many details of this activation are unclear, and have previously been investigated from a computational angle. Yet, these studies have either employed only molecular mechanics (MM), which are inaccurate for metal active sites, or they have described only the active site with quantum mechanics (QM) and neglected the effect of the protein. Here, we employ hybrid QM and MM (QM/MM) methods to investigate the first steps of the LPMO mechanism, which is reduction of of Cu(II) to Cu(I) and formation of a Cu(II)-superoxide complex. In the latter complex, the superoxide can bind either in an equatorial or an axial position. For both steps we obtain structures that are markedly different from previous suggestions, based on vacuum calculations. Our calculations show that the equatorial isomer of the superoxide complex is over 60 kJ/mol more stable than the axial isomer, being stabilized by interactions with a second-coordination-sphere Gln residue, showing a possible role for this residue. Coordination of superoxide in this manner is in agreement with recent experimental suggestions.Comment: 21 pages, 6 figures, 3 tables. 20 pages Supporting Informatio

Genetic analysis shows that Rubus vikensis is a distinct species with a disjunct distribution

Rubus vikensis A. Pedersen ex G. Wendt (sect. Corylifolii) was recently described from a restricted area in north western Scania, Sweden. In this investigation, I show that the same species occurs also on the Onsala peninsula in northern Halland and on a single locality in the middle of Halland. It has 35 chromosomes in all parts of the distribution area. Moreover, I show by random amplified polymorphic DNA (RAPD) analysis that R. vikensis is a distinct and well-defined species, clearly separated from the morphologically similar R. wahlbergii, with which it shares the chromosome number

On the mechanistic origin of damped oscillations in biochemical reaction systems

A generalized reaction scheme for the kinetic interaction of two reactants in a metabolic pathway has been examined in order to establish what minimal mechanistic patterns are required to support a damped oscillatory transient‐state kinetic behaviour of such a two‐component system when operating near a steady state. All potentially oscillating sub‐systems inherent in this scheme are listed and briefly characterized. The list includes several mechanistic patterns that may be frequently encountered in biological system (e.g. involving feedback inhibition, feed‐forward activation, substrate inhibition or product activation), but also draw attention to some hitherto unforeseen mechanisms by which the kinetic interaction of two metabolites may trigger damped oscillations. The results can be used to identify possible sources of oscillations in metabolic pathways without detailed knowledge about the explicit rate equations that apply

Effect of the protein ligand in DMSO reductase studied by computational methods

The DMSO reductase family is the largest and most diverse family of mononuclear molybdenum oxygen-atom-transfer proteins. Their active sites contain a Mo ion coordinated to two molybdopterin ligands, one oxo group in the oxidised state, and one additional, often protein-derived ligand. We have used density-functional theory to evaluate how the fourth ligand (serine, cysteine, selenocysteine, OH−, O2–, SH−, or S2–) affects the geometries, reaction mechanism, reaction energies, and reduction potentials of intermediates in the DMSO reductase reaction. Our results show that there are only small changes in the geometries of the reactant and product states, except from the elongation of the Mo[sbnd]X bond as the ionic radius of X[dbnd]O, S, Se increases. The five ligands with a single negative charge gave an identical two-step reaction mechanism, in which DMSO first binds to the reduced active site, after which the S[sbnd]O bond is cleaved, concomitantly with the transfer of two electrons from Mo in a rate-determining second transition state. The five models gave similar activation energies of 69–85 kJ/mol, with SH− giving the lowest barrier. In contrast, the O2– and S2– ligands gave much higher activation energies (212 and 168 kJ/mol) and differing mechanisms (a more symmetric intermediate for O2– and a one-step reaction without any intermediate for S2–). The high activation energies are caused by a less exothermic reaction energy, 13–25 kJ/mol, and by a more stable reactant state owing to the strong Mo[sbnd]O2– or Mo[sbnd]S2– bonds

Comparison of the active-site design of molybdenum oxo-transfer enzymes by quantum mechanical calculations.

There are three families of mononuclear molybdenum enzymes that catalyze oxygen atom transfer (OAT) reactions, named after a typical example from each family, viz., dimethyl sulfoxide reductase (DMSOR), sulfite oxidase (SO), and xanthine oxidase (XO). These families differ in the construction of their active sites, with two molybdopterin groups in the DMSOR family, two oxy groups in the SO family, and a sulfido group in the XO family. We have employed density functional theory calculations on cluster models of the active sites to understand the selection of molybdenum ligands in the three enzyme families. Our calculations show that the DMSOR active site has a much stronger oxidative power than the other two sites, owing to the extra molybdopterin ligand. However, the active sites do not seem to have been constructed to make the OAT reaction as exergonic as possible, but instead to keep the reaction free energy close to zero (to avoid excessive loss of energy), thereby making the reoxidation (SO and XO) or rereduction of the active sites (DMSOR) after the OAT reaction facile. We also show that active-site models of the three enzyme families can all catalyze the reduction of DMSO and that the DMSOR model does not give the lowest activation barrier. Likewise, all three models can catalyze the oxidation of sulfite, provided that the Coulombic repulsion between the substrate and the enzyme model can be overcome, but for this harder reaction, the SO model gives the lowest activation barrier, although the differences are not large. However, only the XO model can catalyze the oxidation of xanthine, owing to its sulfido ligand

What Is the Structure of the E4 Intermediate in Nitrogenase?

Nitrogenase is the only enzyme that can cleave the strong triple bond inN2. The active site contains a complicated MoFe7S9C cluster. It is believed that itneeds to accept four protons and electrons, forming the E4 state, before it can bind N2.However, there is no consensus on the atomic structure of the E4 state. Experimentalstudies indicate that it should contain two hydride ions bridging two pairs of Fe ions,and it has been suggested that both hydride ions as well as the two protons bind on thesame face of the cluster. On the other hand, density functional theory (DFT) studieshave indicated that it is energetically more favorable with either three hydride ions or with a triply protonated carbide ion, depending on the DFT functional. We have performed a systematic combined quantum mechanical and molecular mechanical (QM/MM) study of possible E4 states with two bridging hydride ions. Our calculations suggest that the most favorable structure has hydride ions bridging the Fe2/6 and Fe3/7 ion pairs. In fact, such structures are 14 kJ/mol more stable than structures with three hydride ions, showing that pure DFT functionals give energetically most favorable structures in agreement with experiments. An important reason for this finding is that we have identified a new type of broken-symmetry state that involves only two Fe ions with minority spin, in contrast to the previously studied states with three Fe ions with minority spin. The energetically best structures have the two hydride ions on different faces of the FeMo cluster, whereas better agreement with ENDOR data is obtained if they are on the same face; such structures are only 6−22 kJ/mol less stable

How are hydrogen bonds modified by metal binding?

We have used density functional theory calculations to investigate how the hydrogen-bond strength is modified when a ligand is bound to a metal using over 60 model systems involving six metals and eight ligands frequently encountered in metalloproteins. We study how the hydrogen-bond geometry and energy vary with the nature of metal, the oxidation state, the coordination number, the ligand involved in the hydrogen bond, other first-sphere ligands, and different hydrogen-bond probe molecules. The results show that, in general, the hydrogen-bond strength is increased for neutral ligands and decreased for negatively charged ligands. The size of the effect is mainly determined by the net charge of the metal complex, and all effects are typically decreased when the model is solvated. In water solution, the hydrogen-bond strength can increase by up to 37 kJ/mol for neutral ligands, and that of negatively charged ligands can increase (for complexes with a negative net charge) or decrease (for positively charged complexes). If the net charge of the complex does not change, there is normally little difference between different metals or different types of complexes. The only exception is observed for sulphur-containing ligands (Met and Cys) and if the ligand is redox-active (e.g. high-valence Fe-O complexes)

Comparison of end-point continuum-solvation methods for the calculation of protein-ligand binding free energies.

We have compared the predictions of ligand-binding affinities from several methods based on end-point molecular dynamics simulations and continuum solvation, i.e. methods related to MM/PBSA (molecular mechanics combined with Poisson-Boltzmann and surface area solvation). Two continuum-solvation models were considered, viz. the Poisson-Boltzmann (PB) and generalised Born (GB) approaches. The non-electrostatic energies were also obtained in two different ways, viz. either from the sum of the bonded, van der Waals, non-polar solvation energies, and entropy terms (as in MM/PBSA), or from the scaled protein-ligand van der Waals interaction energy (as in the linear interaction energy approach, LIE). Three different approaches to calculate electrostatic energies were tested, viz. the sum of electrostatic interaction energies and polar solvation energies, obtained either from a single simulation of the complex or from three independent simulations of the complex, the free protein, and the free ligand, or the linear-response approximation (LRA). Moreover, we investigated the effect of scaling the electrostatic interactions by an effective internal dielectric constant of the protein (ε(int) ). All these methods were tested on the binding of seven biotin analogues to avidin and nine 3-amidinobenzyl-1H-indole-2-carboxamide inhibitors to factor Xa. For avidin, the best results were obtained with a combination of the LIE non-electrostatic energies with the MM+GB electrostatic energies from a single simulation, using ε(int) = 4. For fXa, standard MM/GBSA, based on one simulation and using ε(int) = 4-10 gave the best result. The optimum internal dielectric constant seems to be slightly higher with PB than with GB solvation. Proteins 2012. © 2012 Wiley-Liss, Inc

Will molecular dynamics simulations of proteins ever reach equilibrium?

We show that conformational entropies calculated for five proteins and protein-ligand complexes with dihedral-distribution histogramming, the von Mises approach, or quasi-harmonic analysis do not converge to any useful precision even if molecular dynamics (MD) simulations of 380-500 ns length are employed (the uncertainty is 12-89 kJ mol(-1)). To explain this, we suggest a simple protein model involving dihedrals with effective barriers forming a uniform distribution and show that for such a model, the entropy increases logarithmically with time until all significantly populated dihedral states have been sampled, in agreement with the simulations (during the simulations, 52-70% of the available dihedral phase space has been visited). This is also confirmed by the analysis of the trajectories of a 1 ms simulation of bovine pancreatic trypsin inhibitor (31 kJ mol(-1) difference in the entropy between the first and second part of the simulation). Strictly speaking, this means that it is practically impossible to equilibrate MD simulations of proteins. We discuss the implications of such a lack of strict equilibration of protein MD simulations and show that ligand-binding free energies estimated with the MM/GBSA method (molecular mechanics with generalised Born and surface-area solvation) vary by 3-15 kJ mol(-1) during a 500 ns simulation (the higher estimate is caused by rare conformational changes), although they involve a questionable but well-converged normal-mode entropy estimate, whereas free energies estimated by free-energy perturbation vary by less than 0.6 kJ mol(-1) for the same simulation

Does the DFT Self-Interaction Error Affect Energies Calculated in Proteins with Large QM Systems?

We have examined how the self-interaction error in density-functional theory (DFT) calculations affects energies calculated on large systems (600-1000 atoms) involving several charged groups. We employ 18 different quantum mechanical (QM) methods, including Hartree-Fock, as well as pure, hybrid, and range-separated DFT methods. They are used to calculate reaction and activation energies for three different protein models in vacuum, in a point-charge surrounding, or with a continuum-solvent model. We show that pure DFT functionals give rise to a significant delocalization of the charges in charged groups in the protein, typically by 0.1 e, as evidenced from the Mulliken charges. This has a clear effect on how the surroundings affect calculated reaction and activation energies, indicating that these methods should be avoided for DFT calculations on large systems. Fortunately, methods such as CAM-B3LYP, BHLYP, and M06-2X give results that agree within a few kilojoules per mole, especially when the calculations are performed in a point-charge surrounding. Therefore, we recommend these methods to estimate the effect of the surroundings with large QM systems (but other QM methods may be used to study the intrinsic reaction and activation energies)